Bayesian Modeling for Socio-Environmental Data
Swiss Birds Occupancy Modeling
August 08, 2016
I. Motivation
Modeling presence or absence is a classic problem involving mixture models, specifically random variables that are zero-inflated. Extra zeros are encountered when we model presence or absence because zeros arise from two conditions: truly absent individuals and individuals present but undetected. This means we need a model for the process that controls occupancy, the true state, and model of the data that accounts for detection. This is often our starting point in Bayesian analysis – there a true, unobserved state we seek to understand using a model of a process. We take imperfect observations on that state and must correct them using a model of the data.
II. Problem
A fundamental question in landscape ecology seeks to understand how landscape structure shapes variation in the abundance of species. We will use data from the Swiss Survey of Common Breeding Birds, courtesy of Royle and Dorazio (2008), to model habitat occupancy by a common, resident bird in the Swiss Alps, the willow tit (Parus montanus). The data come from annual surveys of one km2 quadrats distributed across Switzerland (Fig. 1). Surveys are conducted during the breeding season on three separate days, but some quadrats have missing data so that the number of replicate observations is fewer than three.
During each survey, an observer records every visual or acoustic detection of a breeding species (we do not differentiate between these two types of detections in this problem) and marks its location using a global positioning system or, in earlier years, a paper map. Because we are observing a resident species during the breeding season, we assume that the true state (occupied or unoccupied) does not change among sample dates, an assumption known as closure.
Fig. 1. The willow tit (left, c/o of Francis C. Franklin) is one of 70 bird species that are surveyed annually for abundance in 267 1 km2 sampling units distributed across Switzerland (right, c/o the Swiss Ornithological Institute).
We want to understand the influence of forest cover and elevation on the distribution of the willow tit. The data frame SwissBirds has the number of times a quadrat (quadrat) was searched (numberVisits) and the number of times willow tits were detected (numberDetections). We have covariates on forest canopy cover (forestCover) as well as elevation in meters (elevation) for each quadrat surveyed. Develop a model of the influence of forest cover and elevation on the distribution of willow tits. Your model should allow estimation of the optimum elevation of willow tit habitat at the mean forest cover, where optimum elevation is defined as the elevation where probability of occupancy is maximum.
- Diagram the network of knowns and unknowns.
- Write out a mathematical expression for the posterior and the joint distribution.
- Estimate the posterior distributions of the model parameters using JAGS. Some hints:
- You will need to standardize the covariates by subtracting the mean and dividing by the standard deviation for each observation in the elevation and forest cover data. Use the
scalefunction to scale these covariates (this will drastically speed convergence). - You must give initial values of 1 to all unknown 0 or 1 z states.
{ # Extra bracket needed only for R markdown files
sink("SwissBirds.R")
cat("
model{
# priors
p ~ dbeta(1, 1)
for (i in 1:4){
beta[i] ~ dnorm(0, .368)
}
# likelihood
for (i in 1:N){
z[i] ~ dbern(psi[i])
logit(psi[i]) <- beta[1] + beta[2]*elevation[i] + beta[3]*elevation[i]^2 + beta[4]*forestCover[i]
y[i] ~ dbin(z[i] * p, n[i])
}
# derived quantities
elevationMaxScaled <- -beta[2]/(2 * beta[3])
elevationMax <- elevationMaxScaled * sdElevation + muElevation
for (j in 1:length(elevationPredict)){
logit(psiPredict[j]) <- beta[1] + beta[2]*elevationPredict[j] + beta[3]*elevationPredict[j]^2
}
}
", fill = TRUE)
sink()
} # Extra bracket needed only for R markdown fileslibrary(SESYNCBayes)
library(rjags)## Loading required package: coda## Linked to JAGS 4.2.0## Loaded modules: basemod,bugsset.seed(1)
elevation <- scale(SwissBirds$elevation)
muElevation <- attr(elevation, "scaled:center")
sdElevation <- attr(elevation, "scaled:scale")
elevationPredict = (seq(500, 2500, 10) - muElevation)/sdElevation
forestCover <- scale(SwissBirds$forestCover)
muforestCover <- attr(forestCover, "scaled:center")
sdforestCover <- attr(forestCover, "scaled:scale")
inits = list(
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)),
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)),
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)))
data = list(
N = nrow(SwissBirds),
elevation = as.double(elevation),
forestCover = as.double(forestCover),
n = as.double(SwissBirds$numberVisits),
y = as.double(SwissBirds$numberDetections),
muElevation = as.double(muElevation),
sdElevation = as.double(sdElevation),
elevationPredict = as.double(elevationPredict))
n.adapt = 5000
n.update = 10000
n.iter = 10000
jm = jags.model("SwissBirds.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 237
## Unobserved stochastic nodes: 242
## Total graph size: 4128
##
## Initializing modelupdate(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("p", "beta", "elevationMax"), n.iter = n.iter, n.thin = 1)
zj = jags.samples(jm, variable.names = c("elevationMax", "psiPredict"), n.iter = n.iter, n.thin = 1)- Summarize the parameters and check chains for convergence.
summary(zm)##
## Iterations = 15001:25000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] -0.2077 0.27086 0.0015638 0.004115
## beta[2] 1.9649 0.29953 0.0017293 0.004411
## beta[3] -1.0873 0.26438 0.0015264 0.004704
## beta[4] 0.8490 0.22992 0.0013275 0.002841
## elevationMax 1793.0688 144.48064 0.8341594 2.121571
## p 0.7869 0.02871 0.0001658 0.000241
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] -0.7380 -0.3904 -0.2085 -0.02586 0.3294
## beta[2] 1.4046 1.7601 1.9543 2.16092 2.5765
## beta[3] -1.6315 -1.2605 -1.0788 -0.90397 -0.5941
## beta[4] 0.4122 0.6910 0.8441 0.99966 1.3149
## elevationMax 1588.9595 1696.6934 1768.9941 1860.61975 2141.2575
## p 0.7281 0.7680 0.7879 0.80675 0.8399plot(zm)
gelman.diag(zm)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## beta[1] 1 1
## beta[2] 1 1
## beta[3] 1 1
## beta[4] 1 1
## elevationMax 1 1
## p 1 1
##
## Multivariate psrf
##
## 1heidel.diag(zm)## [[1]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.829
## beta[2] passed 1 0.480
## beta[3] passed 1 0.695
## beta[4] passed 1 0.501
## elevationMax passed 1 0.923
## p passed 1 0.384
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.207 0.013792
## beta[2] passed 1.967 0.016177
## beta[3] passed -1.087 0.016369
## beta[4] passed 0.851 0.010028
## elevationMax passed 1794.168 7.289525
## p passed 0.787 0.000811
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.4961
## beta[2] passed 1 0.4231
## beta[3] passed 1 0.8267
## beta[4] passed 1 0.1242
## elevationMax passed 1 0.4175
## p passed 1 0.0789
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.211 0.013468
## beta[2] passed 1.963 0.014393
## beta[3] passed -1.083 0.015020
## beta[4] passed 0.851 0.009627
## elevationMax passed 1793.833 6.962665
## p passed 0.787 0.000805
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.1971
## beta[2] passed 1001 0.4993
## beta[3] passed 1 0.0645
## beta[4] passed 1 0.2333
## elevationMax passed 1 0.1816
## p passed 1 0.1173
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.205 0.014625
## beta[2] passed 1.958 0.014849
## beta[3] passed -1.092 0.016479
## beta[4] passed 0.845 0.009259
## elevationMax passed 1791.205 7.348855
## p passed 0.787 0.000838- What can you conclude about the relative importance of elevation and forest cover in controlling the bird’s distribution? Plot the probability of occupancy as function of elevation at the mean of forest cover. Estimate the optimum elevation of willow tit habitat at the mean forest cover, where optimum elevation is defined as the elevation where probability of occupancy is maximum. Plot a normalized histogram of MCMC output for the optimum elevation at the average forest cover. Overlay 95% credible intervals on the optimum elevation.
par(mfrow = c(1, 2))
psi <- summary(zj$psiPredict, quantile, c(.025, .5, .975))$stat
plot(seq(500, 2500, 10), psi[2,], type = "l", ylim = c(0,1), xlab = "Elevation (m)", lwd = 2,
ylab ="Probability of occupancy")
lines(seq(500, 2500, 10), psi[1,], type = "l", lty = "dashed", lwd = 2)
lines(seq(500, 2500, 10), psi[3,], type = "l", lty = "dashed", lwd = 2)
abline( v = mean(zj$elevationMax), lwd = 2)
text(1300, .9, "Optimum elevation", cex = .8)
hist(zj$elevationMax, breaks = 100, main = "", xlab = "Elevation (m)", xlim = c(1500, 2500),
freq = FALSE, col = "gray")
abline(v = summary(zj$elevationMax, quantile, c(.975))$stat, lty = "dashed", lwd = 2)
abline(v = summary(zj$elevationMax, quantile, c(.025))$stat, lty = "dashed", lwd = 2)
Fig.2. Probability of occupancy at the mean forest cover as a function of elevation (left panel) and the posterior density of the optimum elevation at the mean forest cover (right panel). Dashed give are .025 and .975 quantiles, also known as .95 equal-tailed Bayesian credible intervals.
III. References
Royle, J.A., and R.M. Dorazio. 2008. Hierarchical Modeling and Inference in Ecology: The Analysis of Data from Populations, Metapopulations, and Communities. Academic Press, London, United Kingdom.